Crop Productivity Analysis#

Data#

The following datasets are utilized in this analysis for calculating and mapping crop productivity over the past years:

  1. Dynamic World Dataset:

    • Source: Dynamic World - Google and the World Resources Institute (WRI)

    • Description: The Dynamic World dataset provides a near real-time, high-resolution (10-meter) global land cover classification. It is derived from Sentinel-2 imagery and utilizes machine learning models to classify land cover into nine distinct classes, including water, trees, grass, crops, built areas, bare ground, shrubs, flooded vegetation, and snow/ice. The dataset offers data with minimal latency, enabling near-immediate analysis and decision-making.

    • Spatial Resolution: 10 meters.

    • Temporal Coverage: Data is available since mid-2015, updated continuously as Sentinel-2 imagery becomes available capturing near Real-time.

  2. MODIS Dataset:

    • Source: NASA’s Moderate Resolution Imaging Spectroradiometer MODIS on Terra and Aqua satellites.

    • Description: The MODIS dataset provides a wide range of data products, including land surface temperature, vegetation indices, and land cover classifications. It is widely used for monitoring and modeling land surface processes.

    • Spatial Resolution: 250 meters.

    • Temporal Coverage: Data is available from 2000 to the present, with daily to 16-day composite products.

  3. GDP and Agricultural Data:

    • Source: Official government statistics and international databases.

    • Description: This dataset includes regional GDP figures, agricultural output, export data, and crop production statistics. It is used to analyze the economic impact of agricultural productivity.

    • Spatial Resolution: Varies by dataset; typically available at national and subnational levels.

    • Temporal Coverage: Varies by dataset; typically available annually or quarterly.

  4. Administrative Boundaries:

    • Source: MIMU (Myanmar Information Management Unit)

    • Description: This dataset provides the administrative boundaries for Myanmar at various levels (national, state/region, district, township). It is used for spatial aggregation and analysis of crop productivity and economic data.

    • Spatial Resolution: Varies by administrative level.

Assumptions and Methodological Notes#

This analysis is based on the following assumptions and methodological considerations:

Temporal Definitions and Aggregations#

  1. Fiscal Year Definition: The first quarter (Q1) of each fiscal year starts in April and ends in March of the following year. This aligns with Myanmar’s fiscal calendar.

  2. Quarterly Aggregations: Monthly data are aggregated into quarters (Q1: April-June, Q2: July-September, Q3: October-December, Q4: January-March).

Methodological Assumptions#

  1. Spatial Aggregation: Regional GDP values are based on single-year reference data and are assumed to maintain consistent spatial distributions across all years in the analysis period.

  2. EVI Processing:

    • EVI values are median-aggregated within administrative boundaries to reduce noise

    • Time series preprocessing follows TIMESAT methodology: outlier removal, interpolation, and Savitzky-Golay smoothing

    • Seasonality parameters (start of season, middle of season, end of season) are extracted using a 20% threshold of maximum EVI

  3. Cropland Identification: Pixels are classified as cropland based on Dynamic World’s machine learning classification. Classification accuracy varies by region and land cover type, with potential confusion between cropland and similar land covers.

def find_project_root(marker_file: str = "pyproject.toml"):
    for parent in Path().cwd().parents:
        if (parent / marker_file).exists():
            return parent
    raise FileNotFoundError(f"No project root found with marker file: {marker_file}")


PROJECT_ROOT = Path().cwd().parent.parent
DATA_PATH = PROJECT_ROOT / "data"
EVI_PATH = DATA_PATH / "EVI and Crop Land" / "EVI 2010-2025"
CROPLAND_PATH = DATA_PATH / "EVI and Crop Land" / "Crop Land"
BOUNDARIES_PATH = DATA_PATH / "Shapefiles"

Crop area statistics#

Name PCODE Crop Area (ha) in 2015 Crop Area (ha) in 2025 % Change in Crop Area (2015-2025) Absolute % Change in Crop Area (2015-2025)
0 Sagaing MMR005 897012.27 1355659.42 51.13 51.13
1 Ayeyarwady MMR017 477942.32 1185667.23 148.08 148.08
2 Mandalay MMR010 802809.95 919394.40 14.52 14.52
3 Magway MMR009 739612.49 832402.33 12.55 12.55
4 Bago (East) MMR007 3294.37 612430.37 18490.21 18490.21
5 Yangon MMR013 1962.76 469040.30 23796.98 23796.98
6 Bago (West) MMR008 309351.34 442237.44 42.96 42.96
7 Shan (South) MMR014 4505.87 300314.13 6564.95 6564.95
8 Mon MMR011 7608.89 178099.30 2240.67 2240.67
9 Kachin MMR001 99817.44 164590.06 64.89 64.89
10 Shan (North) MMR015 34504.25 154754.92 348.51 348.51
11 Rakhine MMR012 178160.82 142033.64 -20.28 20.28
12 Nay Pyi Taw MMR018 85037.49 117129.11 37.74 37.74
13 Kayin MMR003 11609.83 82406.40 609.80 609.80
14 Shan (East) MMR016 19354.78 36171.34 86.89 86.89
15 Tanintharyi MMR006 11120.88 12136.31 9.13 9.13
16 Kayah MMR002 0.00 10242.61 0.00 0.00
17 Chin MMR004 6177.81 10230.50 65.60 65.60

Crop Seasonality#

Using this time series dataset of EVI images, we apply several pre-processing steps to extract critical phenological parameters: start of season (SOS), middle of season (MOS), end of season (EOS), length of season (LOS), etc. This workflow is heavily inspired by the TIMESAT software.

Pre-processing steps

  1. Remove outliers from dataset on per-pixel basis using median method: outlier if median from a moving window < or > standard deviation of time-series times 2.

  2. Interpolate missing values linearly

  3. Smooth data on per-pixel basis (using Savitsky Golay filter, window length of 3, and polyorder of 1)

Phenology Process
We then extract crop seasonality metrics using the seasonal amplitude method from the phenolopy package.

The chart below shows the result of this process for a single crop pixel. The blue dots represent the raw EVI values, the black line represents the processed EVI values, and the dotted lines represent season parameters extracted for that pixel: start of season, peak of season, and end of season.

Based on the phenology process, we identified the seasonality to start in July and end in February with the peak being in October. This can vary with geographic region and crop type as well, however, that has not been taken into consideration in this version.

The following figure shows the median EVI in Myanmar from 2010 to 2025 during the crop growing season (June to February). The shaded area represents the interquartile range (IQR) of EVI values across all pixels, indicating the variability in vegetation health during this period.

Harvested Area and Satellite-Derived Crop Area Ranking (2021-2023)#

../../_images/090a422cf226578c6478344c024df1576252fef2453c5915b1988ea570bc284e.png
../../_images/479d0079b02b22d1591242af06f8157ad7c15f8f13c86b32f5a9f22cc0cd811e.png

EVI and Crop Exports#

../../_images/664ae63ab0400c049dc48b227ab7918076c0133b342998ad86466ce1b94197ef.png
../../_images/b449774b77b7694cc5521b8f09fa3563457a8b98423db5dc36ac59e939afc740.png
Dependent variable: total_vegetable_exports_log
(1)(2)(3)(4)
Intercept5.524***5.961***7.502***7.523***
(0.156)(0.208)(0.504)(0.499)
evi_lag_3_log0.348***0.323***0.449***
(0.110)(0.104)(0.120)
evi_lag_6_log0.658***0.408*
(0.192)(0.227)
evi_log-0.448***-0.480***0.0280.064
(0.114)(0.109)(0.189)(0.188)
is_crop_season[T.True]-0.208**
(0.103)
Observations168165162162
R20.0850.1520.2290.249
Adjusted R20.0790.1410.2150.230
Residual Std. Error0.389 (df=166)0.369 (df=162)0.346 (df=158)0.343 (df=157)
F Statistic15.376*** (df=1; 166)14.510*** (df=2; 162)15.669*** (df=3; 158)13.011*** (df=4; 157)
Note:*p<0.1; **p<0.05; ***p<0.01
Total VegBeans/PulsesCornRice
(1)(2)(3)(4)
Intercept7.523***5.800***5.012***6.662***
(0.499)(0.585)(1.469)(0.770)
evi_log0.064-0.107-0.2490.324
(0.188)(0.220)(0.553)(0.290)
evi_lag_3_log0.449***-0.2260.679*0.862***
(0.120)(0.141)(0.354)(0.185)
evi_lag_6_log0.408*0.777***0.4430.232
(0.227)(0.266)(0.668)(0.350)
is_crop_season[T.True]-0.208**-0.080-0.201-0.170
(0.103)(0.120)(0.302)(0.158)
Observations162162162162
R20.2490.3100.0790.147
Adjusted R20.2300.2920.0550.125
Residual Std. Error0.343 (df=157)0.402 (df=157)1.010 (df=157)0.529 (df=157)
F Statistic13.011*** (df=4; 157)17.618*** (df=4; 157)3.351** (df=4; 157)6.751*** (df=4; 157)
Note:*p<0.1; **p<0.05; ***p<0.01
../../_images/1861a0c92dd049ec198f7ce31b98190241f44deee33cb7a43d25e584ba3d1784.png
../../_images/3ea54656345de910ece886a9e37f4aadef237a010e0b7e186fcc63be809ec605.png
OLS Regression Results
Dep. Variable: gdp_overall_log R-squared: 0.050
Model: OLS Adj. R-squared: -0.023
Method: Least Squares F-statistic: 0.6895
Date: Thu, 06 Nov 2025 Prob (F-statistic): 0.421
Time: 19:08:33 Log-Likelihood: 5.7806
No. Observations: 15 AIC: -7.561
Df Residuals: 13 BIC: -6.145
Df Model: 1
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept 19.1420 1.265 15.126 0.000 16.408 21.876
evi_log 0.7596 0.915 0.830 0.421 -1.217 2.736
Omnibus: 1.066 Durbin-Watson: 0.232
Prob(Omnibus): 0.587 Jarque-Bera (JB): 0.678
Skew: -0.495 Prob(JB): 0.712
Kurtosis: 2.677 Cond. No. 58.4


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
../../_images/f72867794a835c27e0f6c229226781d606ac746b690eaeacd77b9b63cf5ae1da.png
../../_images/1be84db0810228f139db721592fb3bd903ad996fc588d163c4e9af3812fc0a4f.png
Dependent variable: gdp_log
(1)(2)(3)(4)
Intercept16.994***32.030***24.636***21.795***
(0.612)(2.857)(3.954)(1.223)
evi_lag_1_log2.334***2.540***1.115***
(0.299)(0.294)(0.132)
evi_lag_2_log-1.557**
(0.607)
evi_log1.701***1.173***-0.009-0.396***
(0.447)(0.311)(0.547)(0.147)
is_crop_season[T.True]1.484***
(0.086)
ntl_sum_log-0.937***-0.642***-0.540***
(0.220)(0.238)(0.086)
Observations60535353
R20.2000.7050.7410.959
Adjusted R20.1860.6870.7190.956
Residual Std. Error0.749 (df=58)0.466 (df=49)0.441 (df=48)0.175 (df=48)
F Statistic14.463*** (df=1; 58)39.098*** (df=3; 49)34.311*** (df=4; 48)283.827*** (df=4; 48)
Note:*p<0.1; **p<0.05; ***p<0.01
../../_images/cedc1706378d51374083977ed30162bfe11b14472e252fa823926c020a654a19.png
adm1_name Agriculture
0 Ayeyarwady 16.0%
1 Chin 0.5%
2 Kachin 2.1%
3 Kayah 0.5%
4 Kayin 2.4%
5 Magway 14.1%
6 Mandalay 8.3%
7 Mon 4.4%
8 Nay Pyi Taw 1.1%
9 Rakhine 3.3%
10 Sagaing 17.0%
11 Tanintharyi 6.9%
12 Yangon 4.9%
13 Bago (East) 5.3%
14 Bago (West) 5.3%
15 Shan (South) 2.7%
16 Shan (East) 2.7%
17 Shan (North) 2.7%
../../_images/f85264086c404721ea4b679316a4448ff20a66bf7e344816369fbef069c2328b.png
../../_images/48d4280a8b553d6dea41cf2c15bc4ae4bc21b4c505c09f4294b8a963a388bc7e.png
../../_images/8c14cefe53bb38d11a5d0119877ccb8e78b7e794d78e5a4609f0c2c049e0160f.png
../../_images/2a16c19c7bb9d33fc2c1966dd5e9b411369ab3a234e885b8a446016808af111b.png
../../_images/200acb6a05321ce27e8407384c93d5153513b77cf9c16e47af94156ee987a8d9.png
../../_images/c5371d70ebebfab8c46682e20aaaa9a9d79a7086dbbd02c2209d2c1fe54e541a.png
Dependent variable: gdp_log
(1)(2)(3)(4)(5)
Intercept11.544***11.028***11.503***6.824***13.923***
(0.176)(0.466)(0.466)(0.405)(0.174)
evi_lag_1_log0.837***-0.0990.674***
(0.148)(0.120)(0.035)
evi_log0.0960.1250.027-1.055***-0.244***
(0.137)(0.149)(0.148)(0.122)(0.036)
is_crop_season[T.True]2.119***1.602***
(0.084)(0.024)
ntl_sum_log0.0530.097**0.150***-0.175***
(0.045)(0.045)(0.035)(0.016)
Observations1057938937937937
R20.0000.0020.0350.4250.963
Adjusted R2-0.000-0.0000.0320.4230.962
Residual Std. Error1.264 (df=1055)1.270 (df=935)1.250 (df=933)0.965 (df=932)0.248 (df=915)
F Statistic0.495 (df=1; 1055)0.871 (df=2; 935)11.275*** (df=3; 933)172.378*** (df=4; 932)1129.960*** (df=21; 915)
Note:*p<0.1; **p<0.05; ***p<0.01
../../_images/2e0c120189c6071656d3225918ca8b8c236ca1a87b3b5830cb657b37609380f7.png